function dvf_pseudo = create_psuedo_dvf_vec_volt3d_IterateVoltage_vectorized(patient, volt, flow, scan, ii,jj,kk)
    %------------------------------------------------------------------------
    %   This file is part of the
    %   5D-Novel4DCT Toolbox  ("Novel4DCT-Toolbox")
    %   DH Thomas, Ph.D
    %   University of California, Los Angeles
    %   Contact: mailto:dhthomas@mednet.ucla.edu
    %------------------------------------------------------------------------
    % $Author: DHThomas $	$Date: 14-Apr-2014 15-23-24 $	$Revision: 0.1 $
    
    
    tic
    % mkdir('Psuedo_dvfs')
    
    dvf_pseudo=zeros(size(ii,2),size(jj,2),3,size(kk,2));
    tic
    %   cd(model_folder)
    
    volt_2d = permute(repmat(volt, [1 patient.dim(2)]),[2,1]);
    flow_2d = permute(repmat(flow, [1 patient.dim(2)]),[2,1]);
    %     k_2d = permute(repmat([1:patient.dim(3)]', [1 patient.dim(2)]),[2,1]);
    %     j_2d = permute(repmat([1:patient.dim(2)]', [1 patient.dim(3)]),[2,1])';
    %      [k_2d,j_2d]=meshgrid(kk,jj);
    [k_2d,j_2d]=meshgrid(kk,jj);
    E(1,:,:) = ones(size(volt_2d));
    E(2,:,:) = volt_2d;
    E(3,:,:) = flow_2d;
    p_max = 3;
    %     p_volt = polyfit(1:length(volt),volt',15);
    %     p_flow = polyfit(1:length(flow),flow',15);
    %     volt_1000 = spline([1:length(volt)],volt,[1:0.01:length(volt)]);
    %     volt_1000 = polyval(p_volt, [1:0.0001:length(volt)]);
    %     flow_1000 = polyval(p_flow, [1:0.0001:length(volt)]);
    dvf_pseudo_u = zeros(patient.dim);
    dvf_pseudo_v = zeros(patient.dim);
    dvf_pseudo_w = zeros(patient.dim);
    
    if patient.par_toolbox>0
        pctRunOnAll javaaddpath E:\ParforProgMonv2\java
        progressStepSize = 1;
        ppm = ParforProgMon('Creating Psuedo Deformation Vector: ',patient.dim(1),  progressStepSize,800, 120);
        
        parfor slice=1:size(ii,2);
            
            
            constant=load([patient.model_params_folder, sprintf('/constant_%d',slice)]);
            constant_vec_u=constant.data(:,:,1);constant_vec_u=constant_vec_u(:);
            constant_vec_v=constant.data(:,:,2);constant_vec_v=constant_vec_v(:);
            constant_vec_w=constant.data(:,:,3);constant_vec_w=constant_vec_w(:);
            
            alpha=load([patient.model_params_folder,sprintf('/alpha_%d',slice)]);
            alpha_vec_u=alpha.data(:,:,1);alpha_vec_u=alpha_vec_u(:);
            alpha_vec_v=alpha.data(:,:,2);alpha_vec_v=alpha_vec_v(:);
            alpha_vec_w=alpha.data(:,:,3);alpha_vec_w=alpha_vec_w(:);
            
            beta=load([patient.model_params_folder,sprintf('/beta_%d',slice)]);
            beta_vec_u=beta.data(:,:,1);beta_vec_u=beta_vec_u(:);
            beta_vec_v=beta.data(:,:,2);beta_vec_v=beta_vec_v(:);
            beta_vec_w=beta.data(:,:,3);beta_vec_w=beta_vec_w(:);
            
            %     gamma=load(sprintf('gamma_%d',slice));
            %     gamma_vec_u=alpha.data(:,:,1);gamma_vec_u=gamma_vec_u(:);
            %     gamma_vec_v=alpha.data(:,:,2);gamma_vec_v=gamma_vec_v(:);
            %     gamma_vec_w=alpha.data(:,:,3);gamma_vec_w=gamma_vec_w(:);
            
            model_u=[constant_vec_u alpha_vec_u beta_vec_u];% gamma_vec_u];
            model_v=[constant_vec_v alpha_vec_v beta_vec_v];% gamma_vec_v];
            model_w=[constant_vec_w alpha_vec_w beta_vec_w];% gamma_vec_w];
            
            %             E=[ones(length(model_u),1)  volt_2d(:) flow_2d(:)];
            dvf_fit_u=zeros(size(jj,2)*size(kk,2),1);%,size(ii,2),
            dvf_fit_v=zeros(size(jj,2)*size(kk,2),1);%,size(ii,2),
            dvf_fit_w=zeros(size(jj,2)*size(kk,2),1);%,size(ii,2),
            
            %%
            E_iterate = E;
            %             E_diff = zeros([p_max, size(E,2), size(E,3)]);
            % F2 = scatteredInterpolant(j_2d(:),k_2d(:),volt_2d(:));
            % F3 = scatteredInterpolant(j_2d(:),k_2d(:),flow_2d(:));
            
            for p = 1:p_max;
                %     tic
                dvf_fit_w=sum(bsxfun(@times, model_w', E_iterate(:,:)))';
                % toc
                % dvf_fit_w1=sum(bsxfun(@times, model_w, E_iterate)')';
                voxel_w=double(bsxfun(@plus,k_2d(:),dvf_fit_w));
                %                 E_diff(p,:,:) = F2( j_2d, reshape(voxel_w, size(j_2d)) )-squeeze(E_iterate(2,:,:));
                
                aa = round(reshape(voxel_w, size(j_2d)));
                aa(aa<1) =1;aa(aa>size(volt_2d,2)) = size(volt_2d,2);
                E_iterate(2,:,:) = volt(aa);
                E_iterate(3,:,:) = flow(aa);
                
            end
            dvf_fit_u=sum(bsxfun(@times, model_u', E_iterate(:,:)))';
            dvf_fit_v=sum(bsxfun(@times, model_v', E_iterate(:,:)))';
            dvf_fit_w=sum(bsxfun(@times, model_w', E_iterate(:,:)))';
            
            %             toc
            
            
            
            dvf_pseudo_u(slice,jj,kk)=reshape(dvf_fit_u,size(jj,2),size(kk,2));
            dvf_pseudo_v(slice,jj,kk)=reshape(dvf_fit_v,size(jj,2),size(kk,2));
            dvf_pseudo_w(slice,jj,kk)=reshape(dvf_fit_w,size(jj,2),size(kk,2));
            
            
            %          if save_files==1;
            %
            % %         iSaveX( sprintf( 'dvf_pseudo_%d', slice ), dvf_pseudo );
            %     end
            ppm.increment();
            
        end
    else
        
        for slice=1:size(ii,2);
            display(sprintf('Scan %d; Generating DVF: Slice %d of %d', scan, slice, size(ii,2)));
            
            constant=load([patient.model_params_folder, sprintf('/constant_%d',slice)]);
            constant_vec_u=constant.data(:,:,1);constant_vec_u=constant_vec_u(:);
            constant_vec_v=constant.data(:,:,2);constant_vec_v=constant_vec_v(:);
            constant_vec_w=constant.data(:,:,3);constant_vec_w=constant_vec_w(:);
            
            alpha=load([patient.model_params_folder,sprintf('/alpha_%d',slice)]);
            alpha_vec_u=alpha.data(:,:,1);alpha_vec_u=alpha_vec_u(:);
            alpha_vec_v=alpha.data(:,:,2);alpha_vec_v=alpha_vec_v(:);
            alpha_vec_w=alpha.data(:,:,3);alpha_vec_w=alpha_vec_w(:);
            
            beta=load([patient.model_params_folder,sprintf('/beta_%d',slice)]);
            beta_vec_u=beta.data(:,:,1);beta_vec_u=beta_vec_u(:);
            beta_vec_v=beta.data(:,:,2);beta_vec_v=beta_vec_v(:);
            beta_vec_w=beta.data(:,:,3);beta_vec_w=beta_vec_w(:);
            
            %     gamma=load(sprintf('gamma_%d',slice));
            %     gamma_vec_u=alpha.data(:,:,1);gamma_vec_u=gamma_vec_u(:);
            %     gamma_vec_v=alpha.data(:,:,2);gamma_vec_v=gamma_vec_v(:);
            %     gamma_vec_w=alpha.data(:,:,3);gamma_vec_w=gamma_vec_w(:);
            
            model_u=[constant_vec_u alpha_vec_u beta_vec_u];% gamma_vec_u];
            model_v=[constant_vec_v alpha_vec_v beta_vec_v];% gamma_vec_v];
            model_w=[constant_vec_w alpha_vec_w beta_vec_w];% gamma_vec_w];
            
            E_iterate = E;
            %             E_diff = zeros([p_max, size(E,2), size(E,3)]);
            % F2 = scatteredInterpolant(j_2d(:),k_2d(:),volt_2d(:));
            % F3 = scatteredInterpolant(j_2d(:),k_2d(:),flow_2d(:));
            
            for p = 1:p_max;
                %     tic
                dvf_fit_w=sum(bsxfun(@times, model_w', E_iterate(:,:)))';
                % toc
                % dvf_fit_w1=sum(bsxfun(@times, model_w, E_iterate)')';
                voxel_w=double(bsxfun(@plus,k_2d(:),dvf_fit_w));
                %                 E_diff(p,:,:) = F2( j_2d, reshape(voxel_w, size(j_2d)) )-squeeze(E_iterate(2,:,:));
                
                aa = round(reshape(voxel_w, size(j_2d)));
                aa(aa<1) =1;aa(aa>size(volt_2d,2)) = size(volt_2d,2);
                E_iterate(2,:,:) = volt(aa);
                E_iterate(3,:,:) = flow(aa);
                
            end
            dvf_fit_u=sum(bsxfun(@times, model_u', E_iterate(:,:)))';
            dvf_fit_v=sum(bsxfun(@times, model_v', E_iterate(:,:)))';
            dvf_fit_w=sum(bsxfun(@times, model_w', E_iterate(:,:)))';
            
            %     cd('Psuedo_dvfs');
            
            %     dvf_pseudo_u1(slice,jj,kk)=reshape(dvf_fit_u1,size(jj,2),size(kk,2));
            %     dvf_pseudo_v1(slice,jj,kk)=reshape(dvf_fit_v1,size(jj,2),size(kk,2));
            %     dvf_pseudo_w1(slice,jj,kk)=reshape(dvf_fit_w1,size(jj,2),size(kk,2));
            
            dvf_pseudo_u(slice,jj,kk)=reshape(dvf_fit_u,size(jj,2),size(kk,2));
            dvf_pseudo_v(slice,jj,kk)=reshape(dvf_fit_v,size(jj,2),size(kk,2));
            dvf_pseudo_w(slice,jj,kk)=reshape(dvf_fit_w,size(jj,2),size(kk,2));
            
            
            %          if save_files==1;
            %
            % %         iSaveX( sprintf( 'dvf_pseudo_%d', slice ), dvf_pseudo );
            %     end
            %     ppm.increment();
            
        end
    end
    
    toc
    dvf_pseudo(:,:,1,:)=dvf_pseudo_u;
    dvf_pseudo(:,:,2,:)=dvf_pseudo_v;
    dvf_pseudo(:,:,3,:)=dvf_pseudo_w;
    
    
    
    
    
    
